/* 
 * File:   main.c
 * Author: Reinhard Weyand
 *
 * Implementation of the fast evaluation algorithm for the coefficients of the cover polynomial
 * as presented in "Computing The Tutte Polynomial In Vertex-Exponential Time" by Björklund et. al. (Appendix p.3-4).
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <unistd.h>
#include "tools.h"
#include <limits.h>
#include <sys/resource.h>

void printCoefficients() { //print coefficients in rectengular format
    int i = 0, j = 0;
    printf("\n ------------------------------------------------------------- \n coefficients in rectengular format:\n(row -> exponent of y, column -> exponent of falling powers of x) \n\n\n");
    printf("   ");
    for (i = 0; i <= n; i++)
        printf("%d ", i);

    printf("\n");
    printf("  "
            " ");
    for (i = 0; i <= 2 * n; i++)
        printf("-");
    printf("\n");

    for (i = 0; i <= n; i++) {
        for (j = 0; j <= n; j++) {

            if (j != 0) printf(" ");
            else printf("%d |", i);
            printf("%d", C_D[i][j]);
            if (j == n)printf("\n");
        }
    }
    printf("\n");
}

int w(int S, int s, int t, int l) {
    int i, j; //positions in submatrix
    int size = getNumOfNodes(S);
    i = getPos(s, S);
    j = getPos(t, S);
    if (!(ELEM(s, S) && ELEM(t, S)))
        return 0;
    if ((s == t) && (l == 0))
        return 1;
    int* subMatrix = (int*) malloc(size * size * sizeof (int));
    int* result = (int*) malloc(size * size * sizeof (int));
    getSubMatrix(S, subMatrix);
    matrixPow(subMatrix, l, result, size);
    int res = matrix(result, i, j, size);
    //   printf("freeing matrices of size %d x %d \n computed result: %d \n",size,size,res);	
    free(result);
    //	printf("freeing submatrix \n");
    free(subMatrix);
    //if(res < 0) printf("ERROR!!\n");
    //	printf("w returns %d \n",res);
    return res;
}

void computeC_X() { //compute number of spanning directed cycles for X according to björklund et.al. by fast moebius inversi
    int S, j, k, s;
    int min;
    for (S = 0; S < pow2; S++)
        c[S] = 0;
    b_flat = (int*) malloc(pow2 * (n + 1) * sizeof (int));
    int size;
    for (s = 0; s < n; s++)
        for (k = 0; k <= n; k++) {
            //	printf("k	=	%d  \n",k);
            for (S = 0; S < pow2; S++) {
                size = popCnt(S);
                B(S, 0) = w(S, s, s, k);
                for (j = 1; j <= n; j++) {
                    if ((1 << (j - 1)) & S) {
                        B(S, j) = (B(S, j - 1) - B(S^(1 << (j - 1)), j - 1));
                    } else {
                        B(S, j) = B(S, j - 1);
                    }
                }
                if (k == size) {
                    c[S] += B(S, n);

                }
            }
        }
    for (S = 1; S < pow2; S++) {
        if (c[S] % getNumOfNodes(S) != 0) printf("PROBLEM!!!");
	//printf("---------------spanning directed cycles for %d = %d--------------\n",S,c[S]);
        c[S] = c[S] / popCnt(S);
        	
    }
    free(b_flat);
}

void computeP_X() { //compute number of spanning directed paths for X according to björklund et.al. by fast moebius inversion	
    int k, S, j, s, t;
    int size;
    a_flat = (int*) malloc(pow2 * (n + 1) * sizeof (int));
    for (S = 0; S < pow2; S++)
        p[S] = 0;
    for (s = 0; s < n; s++) {
        for (t = 0; t < n; t++) {
            for (k = 0; k <= n; k++) {
                //printf("computing P_X for s=%d,t=%d \n",s,t);  		
                for (S = 0; S < pow2; S++) {
                    size = getNumOfNodes(S);
                    A(S, 0) = w(S, s, t, k);
                    for (j = 1; j <= n; j++) {
                        if ((1 << (j - 1)) & S)
                            A(S, j) = (A(S, j - 1) - A(S^(1 << (j - 1)), j - 1)); //f[j](X) = - f[j-1](X\{j}) + f[j-	 
                        else
                            A(S, j) = A(S, j - 1); //f[j](X) = f[j-1](X)
                    }
                    if (k == (size - 1)) {
                        p[S] += A(S, n);
                        	//printf("---------------spanning directed paths for s=%d, t=%d and S=%d = %d--------------\n",s,t,S,p[S]);
                    }
                }

            }
        }
    }
    for (S = 0; S < pow2; S++)
        if (p[S] < 0) printf("ERROR on comptuing spanning directed paths for S = %d \n", S);
    free(a_flat);
}

void computeP_U_S(set U, set S, polynom* subsetPolynom){
	int s,t,k,i,sign;
	int size = popCnt(U^S);
	int size_S = popCnt(S);
	int intermediate;
	int bin;
	
	polynom poly = *subsetPolynom;

	for(i=0; i < n+3; i++)
		poly[i] = 0;

	for(s=0; s < n; s++){
		for(t=0; t < n; t++){
			for(k = 0; k <= size; k++){
				sign = (k % 2 == 0)? 1 : -1;
				intermediate = w(S,s,t,size_S+k-1);
				bin = binCoeff(size,k);
				poly[size_S+k+1] +=  sign*(intermediate*bin);
				
				if(poly[0] < (size_S+k+1)) 
					poly[0] = size_S+k+1;
			}
		}
	}
				
}

void computeP_U_Space(){
	int k, S, j, s, t, U, i, size;
	set* subsets;
	polynom subsetPolynom;
	P_U = (int**) malloc(pow2 * sizeof (int*));
	for(U=0; U < pow2; U++){
		P_U[U] = (int*) malloc((n + 3) * sizeof (int));
		for(i=0; i < n+3; i++)
			P_U[U][i] = 0;
		size = popCnt(U);
		subsets=(int*)malloc(sizeof(int)*myPow(size));
		getSubsets(subsets,U,size);
		for(i=0; i < myPow(size); i++){
			subsetPolynom = (int*) malloc((n + 3) * sizeof (int));
			computeP_U_S(U,subsets[i],&subsetPolynom);
			polyAdd(P_U[U],subsetPolynom,P_U[U]);

			free(subsetPolynom);
			
		}
	
	}
	
}

void computeC_U_S(set U, set S, polynom* subsetPolynom){
	int s,k,i,sign;
	int size = popCnt(U^S);
	int size_S = popCnt(S);
	int intermediate;
	int bin;
	
	polynom poly = *subsetPolynom;

	for(i=0; i < n+2; i++)
		poly[i] = 0;
	if(S==0) return;
	
			for(k = 0; k <= size; k++){
				for(s=0; s < n; s++){
					sign = (k % 2 == 0)? 1 : -1;
						intermediate = sign*w(S,s,s,size_S+k)*binCoeff(size,k);

						poly[size_S+k+1] +=  intermediate;
				
				if(poly[0] < (size_S+k+1)) 
					poly[0] = size_S+k+1;
			}
		}
	
}

void computeC_U_Space(){
	int k, S, j, s, t, U, i, size;
	set* subsets;
	polynom subsetPolynom;
	C_U = (int**) malloc(pow2 * sizeof (int*));
	for(U=0; U < pow2; U++){
		C_U[U] = (int*) malloc((n + 2) * sizeof (int));
		for(i=0; i < n+2; i++)
			C_U[U][i] = 0;
		size = popCnt(U);
		subsets=(int*)malloc(sizeof(int)*myPow(size));
		getSubsets(subsets,U,size);
		for(i=0; i < myPow(size); i++){
			subsetPolynom = (int*) malloc((n + 2) * sizeof (int));
			computeC_U_S(U,subsets[i],&subsetPolynom);
			polyAdd(C_U[U],subsetPolynom,C_U[U]);
			free(subsetPolynom);
			
		}
		for(i=2; i < n+2; i++){
			C_U[U][i] /= (i-1);
		}
		
		
	}
	
}

void computeC_U_moebius() {
    int S, i, size, j;
    C_U = (int**) malloc(pow2 * sizeof (int*));
    int** c_flat = (int**) malloc(pow2 * (n + 1) * sizeof (int*));
    for (j = 0; j < (n + 1); j++) {
        for (S = 0; S < pow2; S++) {
            C_POLY(S, j) = (int*) malloc((n + 2) * sizeof (int));
            C_POLY(S, j)[0] = 0;
        }
    }
    for (S = 0; S < pow2; S++) {
        C_U[S] = (int*) malloc((n + 2) * sizeof (int));
        size = getNumOfNodes(S);
        for (i = 0; i < n + 2; i++)
            C_POLY(S, 0)[i] = 0;
        C_POLY(S, 0)[size + 1] = c[S];
        if (c[S] != 0)
            C_POLY(S, 0)[0] = size + 1;
        for (j = 1; j <= n; j++) {
            if (ELEM(j - 1, S)) {
                polyAdd(C_POLY(S, j - 1), C_POLY(S^(1 << (j - 1)), j - 1), C_POLY(S, j));
            } else {
                polyCopy(C_POLY(S, j - 1), C_POLY(S, j));
            }
        }
        for (i = 0; i <= C_POLY(S, n)[0]; i++) {
            C_U[S][i] = C_POLY(S, n)[i];
        }

    }
    for (i = 0; i < (n + 1); i++) {
        for (S = 0; S < pow2; S++) {
            free(C_POLY(S, i));
        }
    }
    free(c_flat);
}

void computeP_U_moebius() {
    int S, i, size, j;
    P_U = (int**) malloc(pow2 * sizeof (int*));
    int** c_flat = (int**) malloc(pow2 * (n + 1) * sizeof (int*));
    for (j = 0; j < (n + 1); j++) {
        for (S = 0; S < pow2; S++) {
            C_POLY(S, j) = (int*) malloc((n + 3) * sizeof (int));
            C_POLY(S, j)[0] = 0;
        }
    }
    for (S = 0; S < pow2; S++) {
        P_U[S] = (int*) malloc((n + 3) * sizeof (int));
        size = popCnt(S);
        for (i = 0; i < n + 1; i++)
            C_POLY(S, 0)[i] = 0;
        C_POLY(S, 0)[size + 1] = p[S];
        if (p[S] != 0)
            C_POLY(S, 0)[0] = size + 1;
        for (j = 1; j <= n; j++) {
            //		printf("S = %d, j = %d!\n",S,j);
            if (ELEM(j - 1, S)) {
                polyAdd(C_POLY(S, j - 1), C_POLY(S^(1 << (j - 1)), j - 1), C_POLY(S, j));
            } else {
                polyCopy(C_POLY(S, j - 1), C_POLY(S, j));

            }
        }
        for (i = 0; i <= C_POLY(S, n)[0]; i++) {
            P_U[S][i] = C_POLY(S, n)[i];
        }
/*
         printf("***NEW : Path polynomial for S = %d: ", S);
         			for(i=0; i < n+2; i++){
         				printf("%d, ", P_U[S][i]);
         			}
         			printf("\n");
*/
    }
    for (i = 0; i < (n + 1); i++) {
        for (S = 0; S < pow2; S++) {
            free(C_POLY(S, i));
        }
    }
    free(c_flat);
}

void computeC_D(int i, int j) {
    C_D[i][j] = 0;
    int S, vorz, size, k, l;
    int ifac = fac(i);
    int jfac = fac(j);
    for (S = 0; S < pow2; S++) {
        size = getNumOfNodes(S^(pow2 - 1));
        vorz = size & 1 ? -1 : 1;
        
        polyPow(P_U[S], i, P_U_I);
        polyPow(C_U[S], j, C_U_J);
        polyMul(P_U_I, C_U_J, PC_IJ);
        //	printf("done with processing polynoms\n");
        C_D[i][j] += (PC_IJ[0] < n + 1) ? 0 : vorz * PC_IJ[n + 1];
    }
     //   printf("size = %d, C_D[%d][%d] = : %d	**used PC_IJ(n+1) = %d \n", size,i,j,C_D[i][j],PC_IJ[n+1]);


    
    //printf("Dividing %d by %d\n",C_D[i][j],ifac*jfac);
     C_D[i][j] = (C_D[i][j] / (ifac * jfac));
    //printf("faculty product: %d\n",ifac*jfac);

    //  printf("!!! C_D[%d][%d] = %d !!!\n\n", i, j, C_D[i][j]);

}

int main(int argc, char** argv) {
    if (argc < 2) {
        fprintf(stderr, "no file specified\n");
        exit(EXIT_FAILURE);
    }
    int opt = -1;
    char* filename = (char*) malloc(50 * sizeof (char));
    while ((opt = getopt(argc, argv, "f:")) != -1) {
        switch (opt) {
            case 'f':
                strcpy(filename, optarg);
                parseGraph(filename);
                break;
            default:
                fprintf(stderr, "no file specified \n");
        }
    }
    init();
    pow2 = mypow(n);
    module = INT_MAX;
    c = (int*) malloc(pow2 * sizeof (int));
    p = (int*) malloc(pow2 * sizeof (int));
    int S, i, j;
    //computeP_X();
    //computeC_X();
    //computeP_U_moebius();
    //computeC_U_moebius();
    computeP_U_Space();
   computeC_U_Space();
    
    for (i = 0; i <= n; i++)
        for (j = 0; j <= n; j++) {
            computeC_D(i, j);
        }
    printCoefficients();
    free(filename);
    freeAll();
    return (EXIT_SUCCESS);
}
